package twoD.IsoPfem.EmbeddedCrackModel;

import inf.jlinalg.Array2DMatrix;
import inf.jlinalg.ArrayVector;
import inf.jlinalg.BLAM;
import inf.jlinalg.IMatrix;
import inf.jlinalg.IVector;
import inf.jlinalg.MatrixFormat;

import java.util.ArrayList;

import oneD.fem.algorithm.MyLinalg;
import twoD.IsoPfem.CST;
import twoD.IsoPfem.IGeometry;
import twoD.IsoPfem.RealTable;

public class EmbeddedCrackCST extends CST implements IEmbeddedCrackElement {

	private static final double wt = 2e-4;

	private static final int N2_SEPARATE_POSITIVE = 5;
	private static final int N2_SEPARATE_NEGATIVE = 4;
	private static final int N1_SEPARATE_POSITIVE = 3;
	private static final int N1_SEPARATE_NEGATIVE = 2;
	private static final int N0_SEPARATE_POSITIVE = 1;
	private static final int N0_SEPARATE_NEGATIVE = 0;

	private static final int MAX_ITERATION = 100;

	private static final double INTIAL_GAMMA = 0;

	private static final double GAMMA_THRESHOLD = 1e15;

	/**
	 * this variable declare that the element's crack will be fixed. Basically,
	 * if the crack happen at the element, it won't change the direction.
	 * Subsequent computation will compute crack parameters such as crack
	 * opening parameters and compliance parameter
	 */
	private IMatrix H;
	private IMatrix P;
	private double[] crackConfiguration_;

	private boolean isCrack_ = false;// this variable is to indicate cracking
										// condition
	private boolean fixCrack_ = false;
	private double crackSlope_;
	private boolean crackStart_ = false;// this variable turns on when the crack
										// happen at the elements (to impose the
										// initial f0)
	private boolean firstCrack_ = false;

	private double gamma;
	private IVector e;

	private ArrayList<CrackingState> states_ = new ArrayList<CrackingState>();

	public EmbeddedCrackCST(IGeometry T) {
		super(T);
	}

	public EmbeddedCrackCST(RealTable realTable, IGeometry T) {
		super(realTable, T);
	}

	@Override
	public void setCrack(boolean b) {
		isCrack_ = b;
		gamma = INTIAL_GAMMA;
		e = new ArrayVector(2);
		states_ = new ArrayList<CrackingState>();
	}

	@Override
	public boolean isCrack() {
		return isCrack_;
	}

	@Override
	public void checkFirstTimeCrack(LoadingFunction LF, IMatrix Dhat) {
		// check the condition ttilde = ft until it's satisfied
		if (isCrack_ && !firstCrack_) {
			IVector d = getUe();
			IMatrix A = computeA();
			IVector t0 = new ArrayVector(2);
			BLAM.multiply(1.0, BLAM.NO_TRANSPOSE, A, d, 0.0, t0);
			IMatrix Dhati = (IMatrix) Dhat.clone();
			MyLinalg.inverse2x2(Dhati);
			IVector tmp = (IVector) t0.clone();
			BLAM.multiply(1, BLAM.NO_TRANSPOSE, Dhati, t0, 0.0, tmp);
			double ttilde = BLAM.dot(t0, tmp);
			ttilde *= LF.getDnn();
			ttilde = Math.sqrt(ttilde);
			// System.out.println("t0 = " + MatrixFormat.format(t0));
			// System.out.println("Dhati = " + MatrixFormat.format(Dhati));
			// System.out.println("Dnn = " + Dnn);
			System.out.println("Equivalent traction = " + ttilde);
			double ft = LF.getTractionSeparationLaw().getTensileStrength();

			// IVector sigma = computeSigma();
			// IMatrix P = computeP();
			// BLAM.multiply(1.0, BLAM.TRANSPOSE, P, sigma, 0.0, tmp);
			// System.out.println("sqrt(sigma^TPP^Tsigma) = " +
			// Math.sqrt(BLAM.dot(tmp, tmp)));

			crackStart_ = ((ttilde - ft) >= 0);
			firstCrack_ = true; // after this firstCrack will not be touched
								// again
		}
	}

	@Override
	public boolean isFirstTimeCrack() {
		return crackStart_;
	}

	@Override
	public void setCrackSlope(double a) {
		crackSlope_ = a;
	}

	@Override
	public double getCrackSlope() {
		return crackSlope_;
	}

	@Override
	public void fixCrack() {
		H = computeH();
		P = computeP();
		// System.out.println("H = " + MatrixFormat.format(H));
		fixCrack_ = true;
	}

	@Override
	public double[] checkCrackPositiveDirection() {
		if (fixCrack_)
			return crackConfiguration_;
		else {
			crackConfiguration_ = new double[3];
			double[] S1 = computeSigma1();
			double a0 = S1[1];
			double b0 = S1[2];
			// System.out.println("a0 = " + a0);
			// System.out.println("b0 = " + b0);
			IVector xe = this.getXe();
			// compute the coordinate of element center of gravity
			double x0 = (xe.get(0) + xe.get(2) + xe.get(4)) / 3;
			double y0 = (xe.get(1) + xe.get(3) + xe.get(5)) / 3;
			// find the implicit equation of crack line under the form
			// ax+by+c=0; actually the crack line is perpendicular to the
			// principal stress so the principal direction was the normal of the
			// crack line; Presumably the crack line went through the center of
			// gravity of the CST element
			double a = a0;
			double b = b0;
			double c = -a0 * x0 - b0 * y0;
			// checking where the crack line cut through the element. We assume
			// the positive direction is from boundary of crack line point to
			// the side contains one node; to check it we use dot product of two
			// vectors
			if ((a * xe.get(0) + b * xe.get(1) + c)
					* (a * xe.get(2) + b * xe.get(3) + c) > 0) {
				// node 0 and 1 at the same side
				// compute vector cog-node 2
				double v1 = xe.get(4) - x0;
				double v2 = xe.get(5) - y0;
				// check if cog-2 vector and (a0,b0) is in the same side w.r.t
				// the line
				if (v1 * a0 + v2 * b0 < 0)
					// not at the same side then the direction of (a0,b0) is
					// reverse
					crackConfiguration_[2] = N2_SEPARATE_NEGATIVE;
				else
					crackConfiguration_[2] = N2_SEPARATE_POSITIVE;
			} else if ((a * xe.get(0) + b * xe.get(1) + c)
					* (a * xe.get(4) + b * xe.get(5) + c) > 0) {
				// node 0 and 2 at the same side
				// compute vector cog-node 1
				double v1 = xe.get(2) - x0;
				double v2 = xe.get(3) - y0;
				// check if cog-1 vector and (a0,b0) is in the same side w.r.t
				// the line
				if (v1 * a0 + v2 * b0 < 0)
					crackConfiguration_[2] = N1_SEPARATE_NEGATIVE;
				else
					crackConfiguration_[2] = N1_SEPARATE_POSITIVE;
			} else {
				// node 1 and 2 at the same side
				// compute vector cog-node 0
				double v1 = xe.get(0) - x0;
				double v2 = xe.get(1) - y0;
				// check if cog-0 vector and (a0,b0) is in the same side w.r.t
				// the line
				if (v1 * a0 + v2 * b0 < 0)
					crackConfiguration_[2] = N0_SEPARATE_NEGATIVE;
				else
					crackConfiguration_[2] = N0_SEPARATE_POSITIVE;
			}
			if ((crackConfiguration_[2] == N0_SEPARATE_NEGATIVE)
					|| (crackConfiguration_[2] == N1_SEPARATE_NEGATIVE)
					|| (crackConfiguration_[2] == N2_SEPARATE_NEGATIVE)) {
				a0 = -a0;
				b0 = -b0;
			}
			crackConfiguration_[0] = a0;
			crackConfiguration_[1] = b0;
			return crackConfiguration_;
		}
	}

	public IMatrix computeH() {
		if (fixCrack_)
			return this.H;
		else {
			double[] tmp = checkCrackPositiveDirection();
			int i = (int) tmp[2];
			double c = tmp[0];
			double s = tmp[1];
			H = new Array2DMatrix(6, 2);
			int j = i / 2;
			// System.out.println("c = " + c + ", s = " + s + ", i = " + i
			// + ", j = " + j);
			if ((i == N2_SEPARATE_POSITIVE) || (i == N1_SEPARATE_POSITIVE)
					|| (i == N0_SEPARATE_POSITIVE)) {
				H.set(2 * j, 0, c);
				H.set(2 * j, 1, -s);
				H.set(2 * j + 1, 0, s);
				H.set(2 * j + 1, 1, c);
			} else if ((i == N2_SEPARATE_NEGATIVE)
					|| (i == N1_SEPARATE_NEGATIVE)
					|| (i == N0_SEPARATE_NEGATIVE)) {
				// for (int k = 0; k < 3; k++) {
				// if (k != j) {
				// H.set(2 * k, 0, -c);
				// H.set(2 * k, 1, -s);
				// H.set(2 * k + 1, 0, -s);
				// H.set(2 * k + 1, 1, c);
				// }
				// }
				H.set(2 * j, 0, c);
				H.set(2 * j, 1, -s);
				H.set(2 * j + 1, 0, s);
				H.set(2 * j + 1, 1, c);
			}
			return H;
		}
	}

	public IMatrix computeP() {
		if (fixCrack_)
			return this.P;
		else {
			double[] tmp = checkCrackPositiveDirection();
			double c = tmp[0];
			double s = tmp[1];
			IMatrix P = new Array2DMatrix(3, 2);
			P.set(0, 0, c * c);
			P.set(0, 1, -c * s);
			P.set(1, 0, s * s);
			P.set(1, 1, c * s);
			P.set(2, 0, 2 * c * s);
			P.set(2, 1, c * c - s * s);
			return P;
		}
	}

	public IMatrix computeA() {
		IMatrix P = computeP();
		IMatrix C = getCMatrix();
		IMatrix B = getBMatrix();
		IMatrix tmp = new Array2DMatrix(2, 3);
		IMatrix A = new Array2DMatrix(2, 6);
		BLAM.multiply(1.0, BLAM.TRANSPOSE, P, BLAM.NO_TRANSPOSE, C, 0.0, tmp);
		BLAM.multiply(1.0, BLAM.NO_TRANSPOSE, tmp, BLAM.NO_TRANSPOSE, B, 0.0, A);
		return A;
	}

	/**
	 * compute initial gradient f0
	 * 
	 * @param LF
	 * @param Ad
	 * @return f0
	 */
	public IVector computef0(LoadingFunction LF, IVector Ad) {
		IVector f0 = (IVector) Ad.clone();
		double Dnn = LF.getDnn();
		double g0 = LF.getTractionSeparationLaw().g(0);
		MyLinalg.scale(Dnn / (g0 * g0), f0);
		return f0;
	}

	/**
	 * this function return stiffness matrix for unloading case where crack has
	 * been initiated in the element
	 */
	@Override
	public IMatrix computeTangentStiffnessMatrix(IMatrix Dhat) {
		IMatrix ke = super.computeTangentStiffnessMatrix();
		IMatrix H = computeH();
		IMatrix A = computeA();
		IMatrix AH = new Array2DMatrix(2, 2);
		BLAM.multiply(1.0, BLAM.NO_TRANSPOSE, A, BLAM.NO_TRANSPOSE, H, 0.0, AH);
		BLAM.add(1 / gamma, Dhat, AH);
		MyLinalg.inverse2x2(AH);
		IMatrix tmp1 = new Array2DMatrix(6, 2);
		BLAM.multiply(1.0, BLAM.NO_TRANSPOSE, ke, BLAM.NO_TRANSPOSE, H, 0.0,
				tmp1);
		IMatrix tmp2 = new Array2DMatrix(6, 2);
		BLAM.multiply(1.0, BLAM.NO_TRANSPOSE, tmp1, BLAM.NO_TRANSPOSE, AH, 0.0,
				tmp2);
		IMatrix tmp3 = new Array2DMatrix(6, 6);
		BLAM.multiply(1.0, BLAM.NO_TRANSPOSE, tmp2, BLAM.NO_TRANSPOSE, A, 0.0,
				tmp3);
		BLAM.add(-1.0, tmp3, ke);
		return ke;
	}

	/**
	 * this function compute crack parameters when the element is under loading
	 */
	@Override
	public IMatrix computeTangentStiffnessMatrix(LoadingFunction LF,
			IMatrix Dhat) {
		IMatrix ke = super.computeTangentStiffnessMatrix();
		IVector d = getUe();
		IMatrix H = computeH();
		IMatrix A = computeA();
		IMatrix AH = new Array2DMatrix(2, 2);
		BLAM.multiply(1.0, BLAM.NO_TRANSPOSE, A, BLAM.NO_TRANSPOSE, H, 0.0, AH);

		IVector Ad = new ArrayVector(e.getSize());
		BLAM.multiply(1.0, BLAM.NO_TRANSPOSE, A, d, 0.0, Ad);
		IVector t = (IVector) Ad.clone();
		BLAM.multiply(-1.0, BLAM.NO_TRANSPOSE, AH, e, 1.0, t);
		IVector f;
		if (crackStart_) {
			f = computef0(LF, Ad);
			crackStart_ = false;
		} else
			f = LF.dFde(e);

		// System.out.println("ke = " + MatrixFormat.format(ke));
		// System.out.println("H = " + MatrixFormat.format(H));
		// System.out.println("A = " + MatrixFormat.format(A));
		// System.out.println("Dhat = " + MatrixFormat.format(Dhat));
		System.out.println("t = " + MatrixFormat.format(t));
		System.out.println("f = " + MatrixFormat.format(f));
		System.out.println("gamma = " + gamma);

		if (gamma == 0 || BLAM.norm(BLAM.NORM_TWO, e) < wt) {
			// special case
			// gamma = 0: no crack
			// norm(e) < 16e-6: introduce delayed embedded crack
			// do nothing, retain ke
		} else if (gamma == Double.POSITIVE_INFINITY || gamma > 1000) {
			double gamma_star;
			if (gamma == Double.POSITIVE_INFINITY)
				gamma_star = 0.0;
			else
				gamma_star = 1.0 / gamma;
			BLAM.add(gamma_star, Dhat, AH);
			BLAM.multiply(-gamma_star, t, f, AH);

			MyLinalg.inverse2x2(AH);
			IMatrix tmp1 = new Array2DMatrix(6, 2);
			BLAM.multiply(1.0, BLAM.NO_TRANSPOSE, ke, BLAM.NO_TRANSPOSE, H,
					0.0, tmp1);
			IMatrix tmp2 = new Array2DMatrix(6, 2);
			BLAM.multiply(1.0, BLAM.NO_TRANSPOSE, tmp1, BLAM.NO_TRANSPOSE, AH,
					0.0, tmp2);
			IMatrix tmp3 = new Array2DMatrix(6, 6);
			BLAM.multiply(1.0, BLAM.NO_TRANSPOSE, tmp2, BLAM.NO_TRANSPOSE, A,
					0.0, tmp3);
			BLAM.add(-1.0, tmp3, ke);
		} else {
			BLAM.scale(gamma, AH);
			BLAM.add(1.0, Dhat, AH);
			// System.out.println("AH before = " + MatrixFormat.format(AH));
			// System.out.println("t = " + MatrixFormat.format(t));
			// System.out.println("f = " + MatrixFormat.format(f));
			BLAM.multiply(-1.0, t, f, AH);
			// System.out.println("AH after = " + MatrixFormat.format(AH));

			MyLinalg.inverse2x2(AH);
			IMatrix tmp1 = new Array2DMatrix(6, 2);
			BLAM.multiply(gamma, BLAM.NO_TRANSPOSE, ke, BLAM.NO_TRANSPOSE, H,
					0.0, tmp1);
			IMatrix tmp2 = new Array2DMatrix(6, 2);
			BLAM.multiply(1.0, BLAM.NO_TRANSPOSE, tmp1, BLAM.NO_TRANSPOSE, AH,
					0.0, tmp2);
			IMatrix tmp3 = new Array2DMatrix(6, 6);
			BLAM.multiply(1.0, BLAM.NO_TRANSPOSE, tmp2, BLAM.NO_TRANSPOSE, A,
					0.0, tmp3);
			BLAM.add(-1.0, tmp3, ke);
		}
		return ke;
	}

	/**
	 * compute the balance equation (Dhat+gamma*AH)e - gamma*Ad = 0
	 * 
	 * @param Dhat
	 * @param AH
	 * @param Ad
	 * @param gamma
	 * @param e
	 * @return
	 */

	public IVector computeCrackBalanceEquation(IMatrix Dhat, IMatrix AH,
			IVector Ad, double gamma, IVector e) {
		IMatrix K = (IMatrix) Dhat.clone();
		BLAM.scale(1 / gamma, K);
		BLAM.add(1.0, AH, K);
		IVector R = (IVector) Ad.clone();
		BLAM.multiply(-1.0, BLAM.NO_TRANSPOSE, K, e, 1.0, R);
		return R;
	}

	/**
	 * update cracking state using simple fixed-point iteration scheme
	 */
	// @Override
	// public void updateCrackingState(LoadingFunction LF, IMatrix Dhat) {
	// IVector d = getUe();
	// double normd = BLAM.norm(BLAM.NORM_TWO, d);
	// double gammabar;
	//
	// if (crackStart_) {
	// /**
	// * for crack inititation, we need a separate process to compute e0
	// * and gamma0. Since the consistency condition was Dnn x
	// * etilde/gamma = g(0)
	// */
	// // trial calculation
	// double g0 = LF.getTractionSeparationLaw().g(0);
	// double etilde = 0.1 * normd;
	// // double Dss = LF.getDss();
	// double Dnn = LF.getDnn();
	// gammabar = Dnn * etilde / g0;
	// } else
	// gammabar = gamma;
	//
	// IVector ebar = (IVector) e.clone();
	// IMatrix K;
	// IVector R;
	//
	// IMatrix H = computeH();
	// IMatrix A = computeA();
	// IMatrix AH = new Array2DMatrix(2, 2);
	// BLAM.multiply(1.0, BLAM.NO_TRANSPOSE, A, BLAM.NO_TRANSPOSE, H, 0.0, AH);
	//
	// IVector Ad = new ArrayVector(2);
	// BLAM.multiply(1.0, BLAM.NO_TRANSPOSE, A, d, 0.0, Ad);
	//
	// // System.out.println("H = " + MatrixFormat.format(H));
	// // System.out.println("A = " + MatrixFormat.format(A) + ", d = "
	// // + MatrixFormat.format(d) + ", Ad = " + MatrixFormat.format(Ad));
	// // System.out.println("AH = " + MatrixFormat.format(AH));
	//
	// if (gamma == Double.POSITIVE_INFINITY) { // special case
	// MyLinalg.inverse2x2(AH);
	// // System.out.println("AH = " + MatrixFormat.format(AH));
	// BLAM.multiply(1.0, BLAM.NO_TRANSPOSE, AH, Ad, 0.0, ebar);
	// System.out.println("ebar = " + MatrixFormat.format(ebar)
	// + ", gammabar = " + gamma);
	// return;
	// }
	//
	// // System.out.println("initial ebar = " + MatrixFormat.format(ebar)
	// // + ", initial gammabar = " + gammabar);
	//
	// double eps = 1.0;
	// // double eps_old = 0.0;
	// // IVector ebar_old = new ArrayVector(ebar.getSize());
	// // double gammabar_old = 0.0;
	// boolean loop_continue = true;
	// int loop_cnt = 0;
	// do {
	// K = (IMatrix) Dhat.clone();
	//
	// BLAM.scale(1 / gammabar, K);
	//
	// BLAM.add(1.0, AH, K);
	//
	// R = (IVector) Ad.clone();
	//
	// // System.out.println("K = " + MatrixFormat.format(K) + ", R = "
	// // + MatrixFormat.format(R));
	//
	// // solve
	// MyLinalg.solve2x2(K, R);
	// // System.out.println("R solved = " + MatrixFormat.format(R));
	//
	// // update
	// ebar.set(0, R.get(0));
	// ebar.set(1, R.get(1));
	// gammabar = LF.F(ebar);
	//
	// /**
	// * below is convergence criteria processing to avoid numerical
	// * locking
	// */
	// eps = BLAM.norm(BLAM.NORM_TWO,
	// computeCrackBalanceEquation(Dhat, AH, Ad, gammabar, ebar))
	// / BLAM.norm(BLAM.NORM_TWO, Ad);
	// // IVector tmp = new ArrayVector(2);
	// // BLAM.multiply(1.0, BLAM.NO_TRANSPOSE, K, ebar, 0.0, tmp);
	// // MyLinalg.add(-1.0, Ad, tmp);
	// // eps = BLAM.norm(BLAM.NORM_TWO, tmp) / normd;
	// if ((eps < TOLERANCE_RELATIVE)
	// || (gammabar == Double.POSITIVE_INFINITY))
	// loop_continue = false;
	// else {
	// // System.out.println("epsilon_old=" + epsilon_old);
	// // System.out.println("delta of error:"
	// // + Math.abs(epsilon - epsilon_old));
	// // if (Math.abs(eps - eps_old) < 1e-30) {
	// // MyLinalg.add(-1.0, ebar, ebar_old);
	// // loop_continue = BLAM.norm(BLAM.NORM_TWO, ebar_old)
	// // + Math.pow(gammabar_old - gammabar, 2) > 1e-30;
	// // gammabar_old = gammabar;
	// // ebar_old = (IVector) ebar.clone();
	// // }
	// // eps_old = eps;
	// }
	// loop_cnt++;
	// System.out.println("ebar = " + MatrixFormat.format(ebar)
	// + ", gammabar = " + gammabar + ", epsilon = " + eps);
	// } while (loop_continue && loop_cnt < MAX_ITERATION);
	//
	// // if (gamma == Double.POSITIVE_INFINITY) {
	// // System.out.println("gammabar reach infinity");
	// // MyLinalg.inverse2x2(AH);
	// // System.out.println("AH = " + MatrixFormat.format(AH));
	// // BLAM.multiply(1.0, BLAM.NO_TRANSPOSE, AH, Ad, 0.0, ebar);
	// // }
	//
	// e.set(0, ebar.get(0));
	// e.set(1, ebar.get(1));
	// gamma = gammabar;
	//
	// if (gamma > GAMMA_THRESHOLD)
	// gamma = Double.POSITIVE_INFINITY;
	//
	// if (loop_cnt >= MAX_ITERATION)
	// System.out
	// .println("the iteration to compute crack is not converged within "
	// + MAX_ITERATION + " steps");
	// System.out.println("ebar = " + MatrixFormat.format(ebar)
	// + ", gammabar = " + gammabar + ", epsilon = " + eps
	// + ", nloop = " + loop_cnt);
	//
	// states_.add(new CrackingState(gamma, (IVector) e.clone()));
	// }

	@Override
	public void updateCrackingState(LoadingFunction LF, IMatrix Dhat) {
		IVector d = getUe();
		double normd = BLAM.norm(BLAM.NORM_TWO, d);
		double gammabar;
		double Dnn = LF.getDnn();

		if (gamma == 0) {
			/**
			 * for crack inititation, we need a separate process to compute e0
			 * and gamma0. Since the consistency condition was Dnn x
			 * etilde/gamma = g(0)
			 */
			// trial calculation
			double g0 = LF.getTractionSeparationLaw().g(0);
			double etilde = 0.1 * normd;
			// double Dss = LF.getDss();
			gammabar = Dnn * etilde / g0;
		} else
			gammabar = gamma;

		IVector ebar = (IVector) e.clone();
		IMatrix K;
		IVector R;

		IMatrix H = computeH();
		IMatrix A = computeA();
		IMatrix AH = new Array2DMatrix(2, 2);
		BLAM.multiply(1.0, BLAM.NO_TRANSPOSE, A, BLAM.NO_TRANSPOSE, H, 0.0, AH);

		IVector Ad = new ArrayVector(2);
		BLAM.multiply(1.0, BLAM.NO_TRANSPOSE, A, d, 0.0, Ad);

		// System.out.println("H = " + MatrixFormat.format(H));
		// System.out.println("A = " + MatrixFormat.format(A) + ", d = "
		// + MatrixFormat.format(d) + ", Ad = " + MatrixFormat.format(Ad));
		// System.out.println("AH = " + MatrixFormat.format(AH));

		if (gamma == Double.POSITIVE_INFINITY) { // special case
			MyLinalg.inverse2x2(AH);
			// System.out.println("AH = " + MatrixFormat.format(AH));
			BLAM.multiply(1.0, BLAM.NO_TRANSPOSE, AH, Ad, 0.0, ebar);
			System.out.println("ebar = " + MatrixFormat.format(ebar)
					+ ", gammabar = " + gamma);
			return;
		}

		// System.out.println("initial ebar = " + MatrixFormat.format(ebar)
		// + ", initial gammabar = " + gammabar);

		double eps = 1.0;
		// double eps_old = 0.0;
		// IVector ebar_old = new ArrayVector(ebar.getSize());
		// double gammabar_old = 0.0;
		boolean loop_continue = true;
		int loop_cnt = 0;
		do {
			K = (IMatrix) Dhat.clone();

			BLAM.scale(1 / gammabar, K);

			BLAM.add(1.0, AH, K);

			R = (IVector) Ad.clone();

			// System.out.println("K = " + MatrixFormat.format(K) + ", R = "
			// + MatrixFormat.format(R));

			// solve
			MyLinalg.solve2x2(K, R);
			// System.out.println("R solved = " + MatrixFormat.format(R));

			// update
			ebar.set(0, R.get(0));
			ebar.set(1, R.get(1));
			gammabar = LF.F(ebar);

			/**
			 * below is convergence criteria processing to avoid numerical
			 * locking
			 */
			eps = BLAM.norm(BLAM.NORM_TWO,
					computeCrackBalanceEquation(Dhat, AH, Ad, gammabar, ebar))
					/ BLAM.norm(BLAM.NORM_TWO, Ad);
			// IVector tmp = new ArrayVector(2);
			// BLAM.multiply(1.0, BLAM.NO_TRANSPOSE, K, ebar, 0.0, tmp);
			// MyLinalg.add(-1.0, Ad, tmp);
			// eps = BLAM.norm(BLAM.NORM_TWO, tmp) / normd;
			if ((eps < TOLERANCE_RELATIVE)
					|| (gammabar == Double.POSITIVE_INFINITY))
				loop_continue = false;
			else {
				// System.out.println("epsilon_old=" + epsilon_old);
				// System.out.println("delta of error:"
				// + Math.abs(epsilon - epsilon_old));
				// if (Math.abs(eps - eps_old) < 1e-30) {
				// MyLinalg.add(-1.0, ebar, ebar_old);
				// loop_continue = BLAM.norm(BLAM.NORM_TWO, ebar_old)
				// + Math.pow(gammabar_old - gammabar, 2) > 1e-30;
				// gammabar_old = gammabar;
				// ebar_old = (IVector) ebar.clone();
				// }
				// eps_old = eps;
			}
			loop_cnt++;
			System.out.println("ebar = " + MatrixFormat.format(ebar)
					+ ", gammabar = " + gammabar + ", epsilon = " + eps);
		} while (loop_continue && loop_cnt < MAX_ITERATION);

		// if (gamma == Double.POSITIVE_INFINITY) {
		// System.out.println("gammabar reach infinity");
		// MyLinalg.inverse2x2(AH);
		// System.out.println("AH = " + MatrixFormat.format(AH));
		// BLAM.multiply(1.0, BLAM.NO_TRANSPOSE, AH, Ad, 0.0, ebar);
		// }

		e.set(0, ebar.get(0));
		e.set(1, ebar.get(1));
		gamma = gammabar;

		if (gamma > GAMMA_THRESHOLD)
			gamma = Double.POSITIVE_INFINITY;

		if (loop_cnt >= MAX_ITERATION)
			System.out
					.println("the iteration to compute crack is not converged within "
							+ MAX_ITERATION + " steps");
		System.out.println("element looping = " + loop_cnt);

		states_.add(new CrackingState(gamma, (IVector) e.clone()));
	}

	/*
	 * update using Newton-raphson iteration scheme
	 */
	// @Override
	// public void updateCrackingState(LoadingFunction LF, IMatrix Dhat) {
	// IVector d = getUe();
	// double normd = BLAM.norm(BLAM.NORM_TWO, d);
	// double gammabar;
	//
	// if (crackStart_) {
	// /**
	// * for crack inititation, we need a separate process to compute e0
	// * and gamma0. Since the consistency condition was Dnn x
	// * etilde/gamma = g(0)
	// */
	// // trial calculation
	// double g0 = LF.getTractionSeparationLaw().g(0);
	// double etilde = 0.1 * normd;
	// // double Dss = LF.getDss();
	// double Dnn = LF.getDnn();
	// gammabar = Dnn * etilde / g0;
	// } else
	// gammabar = gamma;
	//
	// IVector ebar = (IVector) e.clone();
	// IMatrix K;
	// IVector R;
	//
	// IMatrix H = computeH();
	// IMatrix A = computeA();
	// IMatrix AH = new Array2DMatrix(2, 2);
	// BLAM.multiply(1.0, BLAM.NO_TRANSPOSE, A, BLAM.NO_TRANSPOSE, H, 0.0, AH);
	//
	// IVector Ad = new ArrayVector(2);
	// BLAM.multiply(1.0, BLAM.NO_TRANSPOSE, A, d, 0.0, Ad);
	//
	// // System.out.println("H = " + MatrixFormat.format(H));
	// // System.out.println("A = " + MatrixFormat.format(A) + ", d = "
	// // + MatrixFormat.format(d) + ", Ad = " + MatrixFormat.format(Ad));
	// System.out.println("AH = " + MatrixFormat.format(AH));
	//
	// if (gamma == Double.POSITIVE_INFINITY) { // special case
	// MyLinalg.inverse2x2(AH);
	// System.out.println("AH = " + MatrixFormat.format(AH));
	// BLAM.multiply(1.0, BLAM.NO_TRANSPOSE, AH, Ad, 0.0, ebar);
	// System.out.println("ebar = " + MatrixFormat.format(ebar)
	// + ", gammabar = " + gamma);
	// return;
	// }
	//
	// // System.out.println("initial ebar = " + MatrixFormat.format(ebar)
	// // + ", initial gammabar = " + gammabar);
	//
	// double eps = 1.0;
	// double eps_old = 0.0;
	// IVector ebar_old = new ArrayVector(ebar.getSize());
	// double gammabar_old = 0.0;
	// boolean loop_continue = true;
	// int loop_cnt = 0;
	// if (ebar.get(0) == 0.0)
	// ebar.set(0, 1e-3);
	// if (ebar.get(1) == 0.0)
	// ebar.set(1, 1e-3);
	// do {
	// K = (IMatrix) Dhat.clone();
	//
	// BLAM.scale(1 / gammabar, K);
	//
	// BLAM.add(1.0, AH, K);
	//
	// IVector De = (IVector) ebar.clone();
	// BLAM.multiply(1.0, BLAM.NO_TRANSPOSE, Dhat, ebar, 0.0, De);
	// IMatrix tmp = (IMatrix) Dhat.clone();
	// BLAM.zero(tmp);
	// System.out.println("gammabar = " + gammabar);
	// BLAM.multiply(-1.0 / (gammabar * gammabar), LF.dFde(ebar), De, tmp);
	// System.out.println("tmp = " + MatrixFormat.format(tmp));
	// BLAM.add(1.0, tmp, K);
	//
	// R = computeCrackBalanceEquation(Dhat, AH, Ad, gammabar, ebar);
	//
	// System.out.println("K = " + MatrixFormat.format(K));
	// System.out.println("R = " + MatrixFormat.format(R));
	//
	// // solve
	// MyLinalg.solve2x2(K, R);
	// System.out.println("R solved = " + MatrixFormat.format(R));
	//
	// // update
	// ebar.set(0, ebar.get(0) - R.get(0));
	// ebar.set(1, ebar.get(1) - R.get(1));
	// gammabar = LF.F(ebar);
	// System.out.println("ebar = " + MatrixFormat.format(ebar)
	// + ", gammabar = " + gammabar + ", epsilon = " + eps);
	//
	// /**
	// * below is convergence criteria processing to avoid numerical
	// * locking
	// */
	// eps = BLAM.norm(BLAM.NORM_TWO, R) / normd;
	// if ((eps < TOLERANCE_RELATIVE)
	// || (gammabar == Double.POSITIVE_INFINITY))
	// loop_continue = false;
	// loop_cnt++;
	// } while (loop_continue && loop_cnt < MAX_ITERATION);
	//
	// // if (gamma == Double.POSITIVE_INFINITY) {
	// // System.out.println("gammabar reach infinity");
	// // MyLinalg.inverse2x2(AH);
	// // System.out.println("AH = " + MatrixFormat.format(AH));
	// // BLAM.multiply(1.0, BLAM.NO_TRANSPOSE, AH, Ad, 0.0, ebar);
	// // }
	//
	// e.set(0, ebar.get(0));
	// e.set(1, ebar.get(1));
	// gamma = gammabar;
	//
	// if (loop_cnt >= MAX_ITERATION)
	// System.out
	// .println("the iteration to compute crack is not converged within "
	// + MAX_ITERATION + " steps");
	// System.out.println("ebar = " + MatrixFormat.format(ebar)
	// + ", gammabar = " + gammabar + ", epsilon = " + eps
	// + ", nloop = " + loop_cnt);
	//
	// states_.add(new CrackingState(gamma, (IVector) e.clone()));
	// }

	@Override
	public ArrayList<CrackingState> getCrackState() {
		return states_;
	}

	@Override
	public IVector computeEpsilon(double... xi) {
		IMatrix B = getBMatrix();
		IVector Ue = this.getUe();
		IVector epsilon = new ArrayVector(3);
		BLAM.multiply(1.0, BLAM.NO_TRANSPOSE, B, Ue, 0.0, epsilon);
		if (fixCrack_) {
			IMatrix H = computeH();
			IMatrix BH = new Array2DMatrix(3, 2);
			BLAM.multiply(1.0, BLAM.NO_TRANSPOSE, B, BLAM.NO_TRANSPOSE, H, 0.0,
					BH);
			IVector e_ = new ArrayVector(3);
			BLAM.multiply(1.0, BLAM.NO_TRANSPOSE, BH, e, 0.0, e_);
			MyLinalg.add(-1.0, e_, epsilon);
		}
		return epsilon;
	}
}
